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Abstract 

A model for synchronization of globally coupled phase oscillators including 
"inertial" effects is analyzed. In such a model, both oscillator frequencies 
and phases evolve in time. Stationary solutions include incoherent (unsyn- 
chronized) and synchronized states of the oscillator population. Assuming a 
Lorentzian distribution of oscillator natural frequencies, <?(f2), both larger in- 
ertia or larger frequency spread stabilize the incoherent solution, thereby mak- 
ing harder to synchronize the population. In the limiting case g(£l) = 6(0,), 
the critical coupling becomes independent of inertia. A richer phenomenol- 
ogy is found for bimodal distributions. For instance, inertial effects may 
destabilize incoherence, giving rise to bifurcating synchronized standing wave 
states. Inertia tends to harden the bifurcation from incoherence to synchro- 
nized states: at zero inertia, this bifurcation is supercritical (soft), but it 
tends to become subcritical (hard) as inertia increases. Nonlinear stability is 
investigated in the limit of high natural frequencies. 

PACS numbers: 05.45.+b, 05.20.-y, 05.40.+j, 64.60.Ht 
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I. INTRODUCTION 



The dynamical behavior of large populations of nonlinearly coupled oscillators may de- 
scribe many phenomena in Physics, Biology and Medicine, [jl]-|3[]. In particular synchro- 
nization of mean-field coupled phase oscillators with different natural frequencies is nicely 
illustrated by Kuramoto's well-known and extesively analyzed model [§],[5|. To describe cer- 
tain biological phenomena, inertial effects should be added to this model. In ||, Ermentrout 
revisited the special problem of self-synchronization in populations of fireflies of a certain 
kind (the Pteroptyx malaccae). Compared to observed behavior, the approach to oscillator 
synchronization as described by the Kuramoto model seems to be too fast. Thus a more ap- 
propriate adaptive frequency model has been proposed in |]-|| , where the natural frequency 
of an oscillator is a new independent variable, which is allowed to vary with time. Thus an 
oscillator is described by its phase and frequency. From the mathematical standpoint, the 
new model is governed by a system of coupled second-order differential equations containing 
inertial terms, in contrast to the system of first-order differential equations governing the 
Kuramoto model. Indeed inertia slows down synchronization and this may result in better 
agreement between theory and experimental measurements. Other possible biological appli- 
cations of Ermentrout-type models include after-effects in alterations of circadian cycles in 
mammalians, cf. || 

A different set of applications for oscillators with inertia include power systems described 
by the swing equations ||, or by Hamilton equations [|10|| . An important technologically 
relevant application is the study of superconducting Josephson junctions arrays |IT| , [12| . 
Here inertial terms describe the effect of nonzero electrical capacitance. Such effect is often 



far from being negligible, and it is absent in the Kuramoto model, [13] 
In this paper, we consider the model equations of Ref . M , 



mujj = —ujj + Qj + K rjv sm(ipN — Oj) + £j(t), (1) 

? = 1 N, 



where 8j, ojj and Qj denote phase, frequency and natural frequency of the jth oscillator, re- 
spectively. The natural frequencies are distributed with probability density g(Q), which may 
have a single maximum (unimodal distribution), or several peaks (multimodal distribution). 
The positive parameters m and K are the "inertia" and the coupling strength, respectively. 
The complex order parameter defined by 



r N e 



1 N 



(2) 



measures phase synchronization: > if the oscillators are synchronized and = 
if not. Finally, £ 3 -'s are independent identically distributed Gaussian white noises, with 
= 0, (£,i(t)£,j( s )) — 2D5ij5(t — s). White noise terms were not included in When 
the inertial terms vanish, m — 0, Eqs. (|TJ) and ([]) are exactly the Kuramoto model. 

Typically, iV is very large, and oscillator synchronization is conveniently analyzed in the 
limiting case of infinitely many oscillators. In this limit, models with mean- field coupling 



2 



are described by an evolution equation for the one-oscillator probability density, p(9, u, Q, t), 
14]. For the present model this equation is || 



dp D d 2 p 
dt m 2 Ouj 2 

~ -ttK -oo + n + Kr sin(^ - 6))p] - , (3) 
m ouj o9 

where the order parameter is now given by 

re^ = / / / e ie p(6,uj,tt,t)g(n)dndujd9. (4) 

JO J —oo J — oo 

Equations (||) and (£|) should be supplemented with appropriate initial and boundary data 
(p is 27r-periodic in 9 and has suitable decay behavior as uo — > ±00) plus the normalization 
condition, 

/ / p(9,u,Q,t)dud9 = 1. (5) 

JO J -oo 

Differentiating J 2?r f*™ p{9, 10, f2, t) du) d9 with respect to time, and then using Eq. (§) itself, 
together with periodicity in 9 and decay in u, we find that the left side of (||) is time 
independent. Normalization to unity of the initial probability density then implies @ for 
the solution of (|3|). 

In this paper, we study oscillator synchronization and transition from incoherence to 
synchronization in the model (H) - (§])• The incoherent solution of (|3]) - (|5p (or simply in- 
coherence) is a stationary solution which is independent of 9. This solution asigns equal 
probability to all angles and has r = (no order), so it corresponds to lack of oscillator 
synchronization. There are synchronized solutions which branch off from incoherence as 
the coupling among oscillators is increased. These bifurcations describe the synchroniza- 
tion transitions, which we have analyzed and compared to the corresponding ones in the 
Kuramoto model. Our main results are that inertia: (i) may stabilize incoherence, making 
it harder to synchronize oscillators, and (ii) it may harden the synchronization transition. 
In the Kuramoto model (m = 0) or with oscillators with identical natural frequencies, the 
synchronization transition is soft (supercritical bifurcation), whereas it may become hard 
(subcritical bifurcation) if the distribution of natural frequencies has a nonzero spread (uni- 
modal Lorentzian distribution) or several peaks (e.g., a discrete bimodal distribution). The 
methods we have used in our analysis are similar to those previously employed in the Ku- 
ramoto model ]|15| -[19|| : linear stability of incoherence, bifurcation analysis, high-frequency 
singular perturbations and numerical solutions. An important difference is that now we 
do not have an explicit functional form for stationary solutions (as it was the case for the 
Kuramoto model). This has led us to use mode-coupling expansions of the solution and 
solving the corresponding mode-coupling equations. Solutions of these equations in close 
form are not always accesible, so that we have introduced some closure assumptions. The 
results of these uncontrolled assumptions have been compared to direct simulations or to 
approximate amplitude equations and found reasonable in the limit of small inertia. 

The rest of the paper is as follows. In Section II, we find the incoherent solution and 
study its linear stability for several natural frequency distributions. Results are compared 
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with those obtained in the massless case JT5Hl9|. It is found that the critical coupling 



needed to destabilize incoherence increases with m for "unimodal" frequency distributions 
of the Lorentzian type. The critical coupling is independent of m when g(Q) = 5(Q). In 
this case the time needed to reach synchronization increases as m increases. If g(Q) = 
[S(Q — Qq) + S(Q + O )]/2 (discrete bimodal distribution), the critical coupling may grow 
or decrease with m depending on the values of Qq. In Section III, we construct other 
stationary solutions by two procedures: an amplitude expansion for solutions branching off 
from incoherence and a general expansion in Hermite polynomials which is appropriately 
truncated. An exact analytical solution is obtained if g(fl) = 5(Q) (cf. ||), while analytical 
approximations for small m and Q are available in the general case. We have observed that 
inertia tends to harden the synchronization transition: in the Kuramoto model (m = 0) 
or with oscillators with identical natural frequencies, the synchronization transition is soft 
(supercritical bifurcation), whereas it becomes hard (subcritical bifurcation) in the cases of 
unimodal Lorentzian or discrete bimodal frequency distributions. In Section IV, we obtain 
approximations to stable time-dependent solutions of Eq. (|3]) in the "high frequency limit" , 
Q — ► oo [[TjJ. There are partially synchronized nonlinearly stable solutions of standing 
wave type, as in the Kuramoto model |I7|J18| ). Finally, numerical results are presented in 
Section V, and compared to the approximate or exact solutions of previous Sections. Two 
Appendices at the end are devoted to technical details. 

II. LINEAR STABILITY OF THE INCOHERENT SOLUTION 

The incoherent solution is a ^-independent stationary solution of (^). Its order parameter 
is r = according to (f|). Then Eqs. (|3]), decay as u — > ±oo and the normalization condition 
(0) yield the incoherent solution: 

H K ' ' 27rV 2tiD v ; 

To analyze its linear stability, let us consider a small disturbance about incoherence, 

p(9,u,Q,t) = p (u,n) + eri(9,u,Q,t) + 0(e 2 ), (7) 

where £< 1. Normalization of p(9,uj,Q,t) then implies 

/ / 7]{6,u,n,t)du;d6 = 0. (8) 

JO J-oo 

We now introduce (0) into ([3]) and (f|) and equate like terms in e. To order e, the result is 

at o0 m oui m 2 oto 2 

r2ir r+oo f+oo 

du> 



rj(<j),u,Q,t) 

m JO J-oo J-oo 

x sin(0 - 9) g(Q)dndud(j). (9) 



We now insert a trial solution 
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(10) 



(which is 27r-periodic in 9) into @, thereby obtaining 

d 2 b n m{uj — VL) db n m (1 — m\ — inmuj) b ri 



dw< 



+ 



D 



du 



+ 



D 



nrnK (i5 n i - %5 n _i) ^ . , . 



where we have defined the scalar product 

hoo r+oo . 



ip(u, £l)vjj(uj, f2) g(£l) dVL du. 



(11) 



(12) 



OO J — oo 



Notice that b_ n = b n and that b n = because of the normalization condition (|j). 

Equation (|TT| ) can be transformed into a nonhomogeneous parabolic cylinder equation 
by the following change of variable: 



b n (u; Q, A) = exp 



AD 



w 



Pn(w] fl, A), 



— (u ~n + 2nDi). 



(13) 
(14) 



Inserting ([13]) and (|Tj) into (|ll|), we obtain 



d 2 (3 n 
dw 2 



+ 



w 



m (A + mfi + n 2 D) 



w-2iVmD)' 



A) <J„,1. 



(15) 



(Recall that rfcu = ^jD/mdw when using the definition of scalar product). Let us assume 
now that n 7^ ±1 and that Q is a fixed real number. Then the right hand side of ( Jl5| ) is zero 
and the resulting equation has the following eigenvalues 



X pn (n) = -— -n 2 D-inQ, p = 0,1,2, 



associated to the eigenfunctions 

2 

/3 p , n (w; f2, Ap, n ) = D p (w) = H p 



( w 
W2, 



(16) 



(17) 



which are independent of n and Q. In this formula, D p (w) and H p (x) are the parabolic 
cylinder function and the Hermite polynomial of index p, respectively pO"| , |2~If . The eigen- 
values A Pi „(fi) of fll~6|) , with n = ±1, ±2, . . ., p = 0, 1, . . . and f2 belonging to the support of 
g(Q), constitute the continuous spectrum of the linear stability problem. In fact, a nonho- 
mogeneous linear problem with a homogeneous part given by (|15D cannot be solved for an 
arbitrary source term if A = A Pi „. Notice that the continuous spectrum lies to the left side 
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of the imaginary axis if D > and n ^ 0. Then the "eigenvalues" (|16D have negative real 
parts (and therefore correspond to stable modes). As we have already said, the neutrally 
stable modes with n = have zero amplitude due to the normalization condition (|8|). 

If n = 1, we can solve flTTD by means of an expansion in eigenf unctions D p (w), p = 
0, 1,2, . . .. To obtain the generalized Fourier coefficients of j3i, we multiply both sides of 
(PI) by D p (w) and integrate over w. As Jf^ D p (x)D n (x)dx = \f2ixp\ 5 pn (orthogonality 



condition, cf. §7.711.1 of J2UJ), the result is 



m 

v poo ^(^—WmD) 2 



f°° p\-2- %s/mU ) D rl Ann 

x £ A>M> (is) 



where 



PoO) 



<9po 



m(w- i2y/mP) g _ i (iu _ i2v ^2 
(271)1,0 

Once we have found (3i, we can calculate the scalar product (1, 2" _iv/ "^) Since 
this scalar product appears as a factor in both sides of the resulting expression, we can 
divide by it, thereby obtaining an eigenvalue equation for A: 

-inKy/D ™ 1 [+<» ( f _ ivW ?) 2 n 
1 = £- / eU 1 D pPodw 

/ + OO / , \ 2 

e -(f-*/JSB) 
-00 



In Appendix A, we show that this equation may be rewritten as 

~ 4 7 rv / ^D p! 



X 



^ .dSl, (21) 



+ \ + + D 



where A p (x) is defined as 

Ap{x) = / + °° e-^- tx f dw. (22) 

The result of evaluating this integral is (cf. Appendix A): 

Ap{x) = i p V2^e^x p (x > 0). (23) 
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Inserting fl23|) in we obtain 



p=0 



A + L> + in + *- 



(24) 



As m — > 0, this equation coincides with that obtained for the Kuramoto model, [Tj|. Equa- 



tion (24) can be rewritten in terms of incomplete gamma functions as follows [21,22 

ID 



where pT],|2| 



,.mD 



K 



/oo 
[mD 7(771 (X + D + in),mD) 
-oo 

-7(1 + (A + D + .n)m, mD)] ^J^^ 
(A + in)me mD ^(m(\ + D + in),mD) 



(mDy x+D+iQ > 



xg(n) dn . 



7(0,2;) 



e-*t a_1 dt. 



(25) 



(26) 



^From now on, we analyze Eq. ( f24|) for special frequency distributions: 



(a) Unimodal frequency distribution, g(Q) = 5(Q). 

In this case we show that if ReA = 0, then ImA = 0. Thus, the eigenvalues that may 
acquire a positive real part are real. Then the critical coupling is obtained by setting A = 0. 
By subtracting from Eq. (|4]) its complex conjugate, we obtain 

= Im(A)/(ImA,m,D), (27) 

where 



f(lm\,m,D) = Yl 



p=0 p\ 



(ImX)* + (d + 



(28) 



Notice than the even function f(lm\,m,D) decreases monotonically with Im(A) > 0. On 
the other hand, / tends to zero, as ImA — > +oo, and 

777 

f^—( y mD)- mD -f( y mD,mD)>0, as ImA -> 0. (29) 

Thus / does not vanish at finite values of ImA, and therefore the only solution of (|2"7j ) is 
ImA = 0. Setting now A = 0, Equation (|24]) yields the critical coupling K = K c , 
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2D 



„m D 



.71=0 



' —m D) 



nl 



and therefore 
K r = 2D. 



(30) 



Figures [I] and ^] show the largest eigenvalue A as a function of m and K. To compute A 
numerically, we used from Eq. fl7|), and (|10D , that the amplitude of the order parameter is 



Ce 



(31) 



close to incoherence. Then, the goal is to simulate the evolution of the system, choosing 
the initial condition sufficiently close to the incoherent solution, and obtain numerically the 
amplitude order parameter r (t) . Fig. ^| shows that different eigenvalue curves (for different 
m) intersect the horizontal axis, A = 0, at the same value of K, as expected from Eq. (|30|). 
(b) Unimodal Lorentzian frequency distribution, g(£l) 
Eq. pi) becomes 



e/ir 



K 
2D ( 



mD 



D 



X + D + e 
n + mD 



+ £ 

n=l 



—m D) 



711 



K 



m(A + D + e) + n 
x (m D)- rn(x+D+£) [(m £))-^(A+iJ+«)-i e 

-7(m(A + D + £),mD) I 



m e 



mD 



-m D 



(32) 



An explicit solution for A cannot in general be found. Thus, we consider several limiting 
cases corresponding to physically interesting parameter choices. In the small noise limit, 
D C 1, we consider the cases (z) m = 0(1) fixed, and (ii) mD = 1. It is remarkable that 
the expansion 



7(a,x) 



e x x a 



x' 



To ( a )n+l 



(33) 



where (a)k = a(a + 1) • • • (a + k — 1), k = 1,2,..., holds in both cases. If x = mD — > 0, 
a = mD + m(A + e) > x, (^) holds as a convergent expansion, [2^]. If x = mD = 1, 
and a = 1 + m(A + e) — > oo (with fixed A and e of order 1), (|33D holds as an asymptotic 
expansion, p3|. Inserting (|33|) in Eq. we obtain 



2D _ x - a 
tf 7 ~ + a 



1 + 



x 



+ 



.r 



a + 1 (a+l)(a + 2) 




(34) 



Similarly to the unimodal case, it is possible to prove that A is always real. To this purpose, 
notice that replacing A + e in eq. (|32| ) with A in eq. ( p4|) (setting Q = 0) we obtain the same 
equation. The critical coupling K = K c is then found by setting A = in (|3~4]). In case (i), 
we have 



K c = 2e(me + l) 



2(2 + 3m e) 
2 + me 



D 



O (D 2 ) , 



(35) 



S 



In the limit of vanishing mass, we recover the result K c = 2{D + e) valid for the Kuramoto 



model ||15|| . Another important limit is e = 0, which reproduces the unimodal distribution. 
We find K c = 2D, independently of mass. Thus the spread in frequency distribution plays 
an important role in synchronizing populations of oscillators affected by inertia. 



In case (ii), Eq. fl34|) yields 

~" 2Da(a + l) 



from which 



A 



e + — f 3 ± Vl + 2Km) 
2m v ' 



(37) 



This quantity is always real and vanishes for K = K c = 2e{me + 3) + — . Note that K c 
grows roughly linearly with m. Thus oscillator synchronization is made harder by increasing 
inertia in the limit of vanishing noise. This behavior is slightly different from that described 
in 0. There numerical simulations seemed to show that incoherence remains stable up to 
a critical coupling, which was independent of m. The singular nature of the limit D — > 
makes the cause of this discrepancy unclear, although we should mention that no stability 
analysis was conducted in 0. In the opposite limit m — > oo, A — > —e, and incoherence is 
always stable. 

The stability diagram in the parameter space (s,K) is shown in Figure |3| for m = 0.2, 
and compared to that of the Kuramoto model (m = 0). This diagram is obtained from Eq. 
(|52]) with A = 0, for fixed D and m. In this figure, we have also plotted the evolution of the 
order parameter amplitude for the parameter values marked in the stability diagram by (1) 
to (4). In all cases, the initial condition is taken sufficiently close to the incoherent solution, 
r = 0. 

(c) Bimodal frequency distribution, g(fl) = — fio) + + ^o)]- 

Eq. ([24] ) becomes 



K e mD ~ (-m D) n (n + mD) [m(A + D) + n] _ 

~2D~ ^ n\ [m(A + D) + n] 2 + m 2 Q 2 ~ ' ' 

In the high frequency limit, Qq — * oo, we can find an analytical formula for A by inserting 
the following asymptotic expansion for the incomplete gamma function in 

7(a,x)~ , a — > oo, (39) 

a 

where a = m(A + D + i Q ), and x = mD. The result is 

K f \ 1 \ 

~ T [x + D + iQo + \ + D-in ) ' ^ 

which yields 



\ = -D + ^ + ^i^iml-K 2 . (41) 
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ReA = gives the same critical coupling as the Kuramoto model, K c = AD, for the same 
bimodal frequency distribution |16| . 



Figure |] shows the stability diagram, which is obtained from Eq. (RSI) with ReA = 0, for 
D = 1 and three different mass values m = 0.1, 1, 6. The stability diagram of the Kuramoto 
model (m = 0) is also depicted. Notice how the curves corresponding to differerent masses 
tend to K — AD as Qo — > oo, as expected. Fig. || displays the evolution of the order 
parameter amplitude in different regions of the stability diagram corresponding to m = 0.8 
and D = 1. 

In conclusion, increasing m, D, e (inertia, noise, and frequency spread) makes it more 
difficult to synchronize the oscillator population via stationary bifurcations from incoherence. 



III. MODE-COUPLING EQUATIONS AND STATIONARY SOLUTIONS 



Inspired by the previous linear stability analysis, we shall expand the distribution func- 
tion using a basis of parabolic cylinder functions (or, equivalently, Hermite polynomials) of 
unit mean square norm 



p(9, u, Q, t) 
In (|42T), we have defined 



2ttD\—* 



711 



(42) 



n=0 



ip n (u) = n\ 




(43) 



so that ip n ipp duj 



n\2 n 



5 np . The functions c n ( 



, Q, t) are 27r-periodic in 9, and we have 



2tt 



co(0,fi,t) d6 



(44) 



as it follows from the normalization of p(9, u, Q, t). 

We shall find a system of mode-coupling equations for the coefficient functions c n . Then 
we shall try to find stationary solutions for different frequency distributions g(tt). This 
is not so easy in the general case, so that we shall follow a standard approach: We shall 
recognize in the system of equations for the coefficient functions a particular stationary 
solution corresponding to incoherence. Then we shall try to find a bifurcation equation 
for other stationary solutions branching off from incoherence. We shall see that even this 
requires different approximation schemes in order to succeed. 
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A. Mode-coupling equations 



Let us insert Eq. ( |4"2D into the Fokker-Planck equation (^). We then obtain the following 
hierarchy of coupled partial differential equations for c n (9,Q,t), 



dci 
~dt 



dco 
dt 

— Ccq C\ - 

m 



V2 X 



! D_dci 
m 89 

m 89 



(45) 
(46) 



- T — = -y/nCcn-i c n - Vn + 1 

at 



m 



D dc n+1 
m 89 



Here we have defined the operator 



V m 



8 fi + K r sin(^ - 6) 
89 D 



In terms of the functions c„, the equation for the order parameter becomes 

r2ir r+oo 



r e 



i if) 



e ie c (9,n,t)g(Q) dQ d9. 



r sm{ip — 9) 



2tT /" + 00 



sin(0 — 6) Cq(4>, fl, t)g(Q) dQ i 



(47) 



(48) 



(49) 
(50) 



B. Incoherence and bifurcating stationary solutions 



The incoherent ^-independent solution, c n = C n (Q), depends only on Q. It can readily 
be obtained from the above hierarchy of equations, by ignoring all derivatives. The result is 

1 1 /mW 2 

c n (n) = —^=(-) n n . (51) 

Inserting this into the expansion in Eq. (|42|), we indeed recover Eq. (^). In particular 
setting Q = 0, this yields the simplest incoherent solution, C n = l/(2ir), corresponding to 
the unimodal frequency distribution. Other interesting stationary solutions are partially 
synchronized distributions, which depend on 9. Notice that eq. (|45"D implies that c\ does 
not depend on 9. 

In the following we analyze stationary solutions bifurcating from incoherence for different 
frequency distributions. 

(a) Discrete unimodal frequency distribution, g(Q) = 5(Q). 

Let us look for stationary solutions with finitely many nonzero coefficients c n = for 
n > N and g(Q) = 5(0,). Eq. (g7D yields Cc N = 0, and therefore, c N = K N e % rcos (^- 9 ) 
(K N =constant). Inserting this result in (|47D for n = N — 1, we find £cjv„i = —\fNc^/m. 
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The corresponding solution is not periodic in 6 unless Kn = 0. Repeating this argument, 
we obtain the stationary solution 



c o (0) 



e § rcos(^-e) 



for n > and the normalization condition has been used. 

2 7 rD\ 1 /4 



where c n 

the resulting distribution, p(#, cj) = 
to its two arguments, 6 and uj, cf 



(52) 



Observe that 



cq(6)iI> (u) 



-muj 2 /4D 



is factorized with respect 



19] 



It is remarkable that Cq(9) is independent of m, 
and coincides with that obtained for the Kuramoto model. Therefore, from (fi9| ) the order 
parameter does not depend on inertia, and the bifurcation diagram for r is exactly the same 
as that for the Kuramoto model. 



(b) General frequency distributions. 

For general frequency distributions, we could try to find stationary solutions of the mode- 
coupling equations (f4^ ) which bifurcate from incoherence. It is however more direct to work 
with the stationary Fokker-Planck equation (|3p as follows. We consider stationary solutions 
as functions of a fixed value of the synchronization parameter r. Then we expand these 
solutions in power series of r and write a hierarchy of equations for the coefficient functions. 
The first coefficient function should be the incoherent solution of synchronization parameter 
r = 0. Inserting the power series probability density function into Eq. (|4]), we find the 
amplitude equation for stationary solutions bifurcating from incoherence. This procedure is 
explained in Appendix B. We quote here the result: 



2D 



6 



(53) 



where 



a 



1 _L ^ 2 
-°° 1 + -JyZ 



g(Q) olQ 



+ g(-mZW + ^) 2 



P =i 



pi 



x 



-oo (1 



P \2 _i_ ft 2 
mD> D 2 



g(Q) olQ 



(54) 



and j3 can be calculated numerically by solving a system of differential equations. As the 
inertia vanishes, mD —>■ 0+, a becomes 



a 



g({l)dtt m 
i , n 2 7) / i i n 2 

1 + ^2 U J -°° 1 + D 2 



g(Q)dn + 0(m 2 D 2 ). 



(55) 



Notice that a in Equation (^) coincides with 2D/K in Eq. ( p4[) provided A = and 
g(Q) = g(—Q). This means that the critical coupling for bifurcation towards stationary 
synchronized states is obtained at Ka = 2D, no matter what the symmetric frequency 
distribution may be. 
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Before interpreting the amplitude equation (p3|) , we shall outline a procedure to obtain 
its coefficients based upon an uncontrolled closure assumption which is accurate for small 
values oimD. Consider the expansion (f42|). We expect that the coefficients c n of stationary 
solutions close to incoherence do not differ much from the coefficients of the latter. In view 
of the functional form (0), we anticipate that the coefficients c n approach zero as n — > oo 
faster for smaller values of the mass. Thus we shall now consider stationary solutions for 
general frequency distributions, such that c n = in (H7H, for all n > 3. By solving equation 



for such a stationary solution with n = 2, we obtain 



m 

c 2 (9, n) = x / — [n + Kr sin(^ - 9)} c x (0). (56) 



We now insert this expression into (|46j ). The resulting equation is solved for a Co, which is 
27r-periodic in 9 and obeys the normalization condition (|5|). We find 

»<'- n >~ e *iy n) < w > 

y>(0, Q) = /^[l - raKr cos(^ - - rf)] 
Jo 

Xe -±[Qr,+Krco S ^-e-r,)] drj ^ 

Z(O) = cos W- 9 V(0, fi) d0, (59) 

JO 

(60) 



In (|57|) , r should be determined so that 

r2n poo 



r= / co(0,n)cos(ip-6)g(Q,)dnd6, (61) 

JO J-oo 

holds. The function Cx(Q) can be obtained by integrating ([IB|) with respect to 6 and using 
the normalization condition for cq together with the 27r-periodicity in 9 of cq and c 2 - 

1 lm r r^n 
c 1 = —J—tt + Kr sm(ip-9)c (9,tt)d9 . (62) 
2n V D l Jo 

Then, the stationary distribution can be approximated by 

\ m J 

+c 2 (9,Q)M^)]e- mv2/iD . (63) 

Notice that a nonvanishing Ci(f2) in ( f42|) implies that the probability density is no longer 
even in u (ipi is an odd function of uo). By fl6"2|) , this occurs in the synchronized phase for 
the case of nonidentical oscillators. For the bimodal frequency distribution, Fig. ^| shows 
Ci as function of VLq for two different values of m. Both the approximate expression (|62"D 
and results of direct numerical simulations are depicted. Notice that the agreement between 
our approximation and the numerical result improves as m decreases. Thus we observe that 
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for each fixed nonzero Q and each fixed 6, the distribution function is no longer peaked at 
uj = 0. However, the instantaneous frequency distribution, defined by 

/•2tt poo 

p (6,lu,0) g(Q)dQ d0, 



may turn out to be even in uj. This is certainly true for the approximate stationary distri- 
bution (|63|), for 

m Kr r 2jT r°° 



/oo ITT) K V f /"Cxj 

c x {n) g {n)dn = J — — / Co (e,n)sm(^-9) g (n)dnd9 = o, 
-oo V D Z7T JO J-oo 



as it follows from the definition of the order parameter and the expansion (42f). 

Another indication that the exact instantaneous frequency distribution may be even in 
uj is that the average frequency tends to zero as t — > +00. This would occur if the stable 
stationary distribution is even in uj (although, admittedly, distributions which are not even in 
uj may have zero mean). The result can be shown directly from the Fokker-Planck equation 
(0). We multiply that equation by uj g(Vt) and integrate with respect to all variables. Then 
we obtain the following equation for the mean value (uj): 



+ M = f°° n g (n)dn. 

dt m J-00 



The right hand side of this expression is zero if g(—£l) = g(ty- Then (uj) tends to zero 
exponentially fast as t — > +00. 

Let us now go back to the problem of obtaining the coefficients in the amplitude equation 
(|53|). Let us expand fl5T|) in powers of r and insert the result in fl6"ip . We obtain an 
approximation to the coefficients of the amplitude equation ( [Si]) , a is again given by ([53D , 
and the expression for (3 is 



o n 2 _ 3 1 39™f^ 

°D 2 2 4 D 



1 + (4+f)^ 3 

g(n)dQ. (64) 



3m»n> (l + 



We now interpret the amplitude equation (^) in the usual way. Notice that the nonzero 
solution is approximately given by 



Kr 



I3(2D-Ka) 



K[3D 



Then the critical value of K is K* = 2D /a, and the sign of (3 determines the direction of the 
bifurcating branch of stationary solutions. Assume a > 0. Then the bifurcating stationary 
solution exists for K > K* if (3 < (supercritical bifurcation) . If (3 > 0, we have a subcritical 
bifurcation and the partially synchronized stationary solution exists for K < K*. The 
critical coupling tends to infinity as a — > 0+, and the bifurcation does not occur at positive 
couplings if a < 0. Given the formulas ( CT ) or Q55D, a could become negative, depending on 
the value of m. Let us now analyze the bifurcation diagram corresponding to two different 
frequency distributions. 
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Unimodal Lorentzian frequency distribution. 

For an unimodal Lorentzian frequency distribution, the coefficients a and (3 of fl55|) 
and (^) are approximately given by 

a = ——(l-me), (65) 
D + 7 

^4 (Z) + ^(2 C + £ ) '- 1 + ^< 30 + ^ (66) 

Apparently, the sign of a could again be negative provided me is sufficiently large. 
However, evaluation of the exact expression (|54]) shows that a > 0. Figure [7] shows a 
as a function of e for two different masses. Notice that the approximate expression for 
a fits better the numerical result as m and e decrease. 

In contrast with a, (3 can really change sign for a Lorentzian frequency distribution. 
Fig. |8| shows the good qualitative agreement between the approximate expression for 
/3 and its numerical evaluation from the exact equations (P15|) . 

There is a critical mass m c (e and D are kept fixed) for which (3 = 0. This mass 
separates the supercritical and subcritical bifurcation regions. Setting (3 = in eq. 
( |66"D yields 

m c = — -. (67) 



Figure || shows the qualitative agreement between the approximate expression 
and the numerical result obtained from the solution of ( |B15| ). 

Note that in the massless case, the stationary solution always branches off supercrit- 
ically, independently of e. In the presence of inertia and for appropriate values of e 



we have found a subcritical bifurcation. This prediction is illustrated by Figures [LO 
(parameters corresponding to a supercritical bifurcation) and [Tl] (parameters corre- 
sponding to a subcritical bifurcation). Figure [TO is a typical supercritical bifurcation 



diagram (order parameter amplitude versus K). In the subcritical case, Figure |IT 
shows the different evolution of the order parameter amplitude for two different initial 
conditions. Notice the bistability between incoherence and the stable synchronized 
solution typical of a subcritical bifurcation. 

Bimodal frequency distribution 

For the bimodal frequency distribution, a may be zero. This means that there is a 
critical frequency ^g(m) above which K* — oo. In fact, we have shown in Fig. |], 
looking to the branch corresponding to A = 0, that in case of a bimodal frequency 
distribution, it is possible to find some values of fl without its corresponding K in 
the stability diagram of incoherence, due to the existence of an horizontal asymptote 
in the space parameter (K, Qq). This means that there is not a finite coupling K* 
where the stationary solution branches off. On the other hand, this never occurs in 
the bimodal Kuramoto model [m = 0, K c j D = 2(1 + fl^/ D 2 )}, because in this case 
there is not such an horizontal asymptote. 
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Eq. fl55|) shows that the critical frequency at which a = is 

^ = A- («) 



Figure [12] compares the previous approximate expression to the exact critical frequency 
obtained from (|54|). Note that the approximation improves as m decreases. Fig. |I3| 
shows that K* increases as Qo does. K* becomes infinity for Q > Q^. 

For a > 0, there is another important critical frequency. In this case, the sign of (3 
decides whether the bifurcation is subcritical (/3 > 0) or supercritical ((3 < 0). The 
sign of j3 depends on m and VL Q . For small mass, we can use the approximations fl55|) 
and (^) (ignoring terms which are quadratic in the mass). Then we find that the 
critical frequency at which (3 = is 



D 



75 p:- (69) 



8 



We have solved numerically the system of equations ( [B15|) , and compared with the 
approximation (F3[) for f3. In Figure IT3I the coefficient j3 is plotted as a function of f2 



for different masses. Similarly, Figure [L^ shows how the critical frequency Q$ varies as 
a function of m. Note that the analytical approximation improves as m goes to zero 
(as expected). 

Notice that Qq < Q^. Then the branch of synchronized stationary states bifurcates 
subcritically from incoherence at K* = oo, provided Qq > Q^. 



Summarizing, inertia favors the subcritical character of bifurcations describing the tran- 
sition from incoherence to the partial synchronized state. In fact in the bimodal case, the 
transition of incoherence to synchronization will most likely occur as a subcritical bifurcation 
as inertia increases. For a continuous unimodal Lorentzian frequency distribution, inertia 
may turn subcritical the supercritical bifurcation to stationary synchronized states, which 
is always found in the masslesss Kuramoto model. 



IV. HIGH FREQUENCY LIMIT 

The high frequency limit of Eq. (|3|) can be analyzed by means of a multiscale method. 
For the Kuramoto model 0, this method lead to the result that the probability density is (to 
leading order) a superposition of different components corresponding to the different peaks 
of the oscillator frequency distribution. To apply this method to the present model, we shall 
assume that the frequency distribution g(Q) has m maxima located at flo^i, I = 1, m, so 
that g(Q)dfl tends to the distribution 

m 

r(i/) = 5>,<y(i/-n)di/, (70) 
i=i 

as Qq — > oo. In order to simplify the calculations, we change variables to a comoving frame: 
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V 



p = e-nt = e--t, 

e 



and let uj' = uj — Q, where 



(71) 



(72) 



Then we obtain the following equations: 

D 



duo 12 mduj' {P ' 3) W 



dt 



Uj = —uj' + Im K I ^2 a i e 



dp 



(73) 



x /r tJ >p((3',uj',t,vi;e)df3'duj' 



OO p2n 



oo JO 



Pj(P', t, J\ e)dp'duj' 



(74) 
(75) 



where pj = p(P, t, uj', Vj\ e), and p ~ Yl, ajPjS(u — Vj), |I9| . We now define fast and slow time 
scales, t = t/e and t, and make the Ansatz: 



n=0 

Inserting this into the governing equations, we obtain the following hierarchy: 

0. 



(76) 



dpf 



dp? 



dp 



2 Jo) 



<9r ^ 



(0) 



d 2 p 



+ 



Id,, „■„. 



(77) 



—a; 



7 {pf [im (a,e-^z] 



(o) 



dp m du 



+ ^a j e i ^- v ^e- if, z} 0) 



(78) 



where 



z?\t) = j_ 



oo r2n 



-oo JO 



(77, uj, t, uj) drj du. 



(79) 



Elimination of secular terms yields the following condition: 



dp 



(°) 92 (0) qM 
3 I) 1 j .' I j 



dt 
]__d 

m du 



UJ 



du' 2 dp 
- [pf [-uj' + Ka 3 Im (e^zf)] } = 0. 



(80) 
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Notice that this equation corresponds to the equation of a unimodal frequency distribution, 
already studied in ||. Its solution evolves towards the following stationary state as time 
evolves: 



D (K aj/D)RjCos(^j-/3) 



where 



\ e iri pf\7],u)d V duj= Urn Zf{t). (82) 

-oo JO t^oo 

Incoherence of a component corresponds to i?j = 0. We know that as Katj surpasses 
2D, a stable synchronized solution bifurcates from incoherence. In the particular case of a 
symmetric bimodal frequency distribution, the bifurcation value is K = 4D, independently 
of the inertia m. This agrees with the results of the linear stability analysis for Q — > oo. 
All the results previously obtained for the Kuramoto model can be applied to the present 
model without any modification ||. Thus, both a stable standing wave solution (SW) and an 
unstable travelling wave solution (TW) bifurcate supercritically from incoherence at K = 4D 
18| , 19] . In Figure [16] illustrates the comparison between the asymptotic solution (|8~ID and the 



result of direct numerical simulation for a large enough value of the frequency CIq. Finally, 



Fig. T7 shows the global bifurcation diagram of the bimodal case for the case of positive a 
and f3. As explained in the previous Section, the branch of synchronized stationary states 
bifurcates subcritically from incoherence at K* = oo, provided Qq > Q™. Thus K* and the 



end of the SW and TW branches in Fig. 17 should extend to infinity as Qq — > +oo. 



V. NUMERICAL RESULTS 

Four different numerical methods have been used. Numerical simulations of the system of 
Langevin equations ([]]) were carried out for a large number of oscillators (N=20 000), using 
an Euler method. A standard finite difference method was used to solve numerically the 
Fokker-Planck equation (|3|), or the system of partial differential equations ([17]). In addition, 
we have used a simple spectral method, which generalizes the one proposed in . The idea 
is to solve a set of ordinary differential equations for moments of p: 

(xi) k := r Ci (e,n k ,t)cos[j&-e)}de, (83) 
Jo 

id)k-= / 27r Q(^,fi fc ,t)sin[j(^-^)]^, (84) 

27T POO POO 

/ / p(0,w,fi,f)cos(V>-0) 

J —oo J —oo 

xg(Q)dQdud6 . (85) 



The coefficient functions Cj(0, f2, t) are the same as in the parabolic cylinder expansion ( [42] 
The integral in fl85|) will be approximated by a suitable quadrature formula, as 
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J2 a <i c o(^ g , 0, t) cos(t/> - 9) d9 = a q( X l)q- 

o=l J ° q=l 



(86) 



For instance, the Gauss-Laguerre quadrature has been chosen case for the case of a 
Lorentzian frequency distribution. The system of ordinary differential equations is given 
by: 



(M)k = jViJ—(Vi-i)k + v 7 ^ J— 
V m VmD 



fifc(^i_i)/t 



+ 1 



' i K r 
mD 2 



(Vi-i)k - (yLi)k - — (xi) k 



j- 1 



rn 



+ IV Dmj{yi +1 ) k - jip{y-)k, 
^k{yLi)k 



{vl)k = -jViJ— Oi-i)& + v 7 ^ J— 
V m VmD 



/ i Kr r 



)fe (^-l)*; 



- + lVDmj (x J l+1 )k+jtp(xl)k, 
i = l,...,N, j = l,...,M, 



(87) 



(88) 



Oo)fc 




- ji>(vo)k, 



(Vo)k = -j\ —(x{)k+ jip(xo)k, 
m 



0, J 



M. 



(89) 
(90) 



In order to numerically simulate this system, it is necessary to truncate the hierarchy 
after a reasonable number of modes N, and M. These numbers will depend on the inertia 
m, the spread, and the coupling strength. They should be chosen large enough, so that the 
numerical results do not depend on N and M. 



Fig. |18| shows the ratio c™ ax /c™ ax as a function of n, for the stationary solution and for 
various mass values. c™ ax is the maximum value of the stationary coefficient c n (9, flo = 1). It 
has been obtained from a finite-difference solution of Eq. (^). Note that this ratio increases 
as inertia does. Thus the truncation approximation ( |53| ) ceases to make sense for larger 
values of inertia. Whether truncating (at higher order) the moment equations ( p3|) and 
(jS4]) yields good numerical results, is tested in Fig. [19]. This figure shows how closely the 
previous system (with moments of order 4 or 10) resembles the direct solution of the Langevin 
equations. We notice that the system containing moments of order 10 is rather close to the 
solution of the Langevin equations, but it does not contain the fluctuations unavoidable in 
stochastic methods. If we are interested in solving the nonlinear Fokker-Planck equation, 
the system of moment equations seems a good alternative to solving Langevin equations for 
a large number of oscillators. 
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VI. CONCLUSIONS 



We have investigated synchronization properties of a model of globally and nonlinearly 
coupled phase oscillators, where the effects of white noise, inertia and spread in the natural 
frequency distribution are all considered. The linear stability of the incoherent solution is 
rigorously analyzed. Stationary and time-dependent solutions of the standing and travelling 
wave type are obtained by a variety of perturbative methods. These include finding and 
approximately solving mode-coupling equations for expansions of the probability distribution 
in parabolic cylinder functions, finding amplitude equations near bifurcation points, and 
hierarchy-closure assumptions. Numerical simulations of the different equations and the 
original model favorably agree with the different perturbative results. 

Inertia changes the stability boundaries of the incoherent solution in a non trivial way. In 
the case of an unimodal Lorentzian frequency distribution, incoherence is stabilized, but the 
effect of mass completely drops out if the frequency spread vanishes. For a discrete bimodal 
frequency distribution, both stability boundaries and the character of the transition from 
incoherence to synchronized states depend on the values of the natural frequency f2 and 
on inertia. The effect of inertia on the stationary solution is dramatic in some cases. In 
general, inertia hardens the synchronization transition: it may render subcritical (hard) 
originally supercritical (soft) transitions (in unimodal Lorentzian frequency distributions), 
or it increases the region in parameter space where the transition is subcritical (discrete 
bimodal frequency distributions). Analytical as well as numerical calculations confirm these 
findings. 
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APPENDIX A: INTEGRALS OF SPECIAL FUNCTIONS AND EIGENVALUE 

EQUATION 

In (p0|) there appear two integrals over w: 
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.4 



f + °° e^- l ^) 2 D pP ' dw, 

J —oo 




The second integral is directly A P {\J mD). The first integral equals 



(2n)?D J-oo 
because of fll~9|) . An equivalent expression is 
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This can be written as 

(2n)iD J-oo PK ' dw 

m r+oo d 



i(2n)2D J- 



( °° DJw)^-e~^ w ~ i2x)2 dw 

J-00 OX 



provided we set x = ymD after performing the calculations. Thus we see that 



i{2-k)iD 

This expression and the previous one for B yield (|2T|). 

Let us now obtain A p (x) for integer p and x > 0. By using that D p (x) is even for even 
p and odd for odd p, we can write 

Ap(x) = e x2 / + °° e ~^ +lwx D p (w) dw 
J — 00 

= e x e 4 (cosxw; + ismxw) D p {w) dw 

J —00 

= 2e x / DJw) e 4 f{xw) dw, 
Jo 

where f(xw) = cosxw for p even and f(xw) = ismxw for p odd. We can find (cf. §7.741 of 

my- ' 



r+00 2 rzf 2 

J D p {w) e 4 cosxw dw = (— 1) y 7^ e 2 x 

/•+OO m 2 ASF 

y D p (w) sinxwdio = (— l)y — 
provided x > 0. From these formulas, we obtain 



2n 



e-Tx 2n+1 



An(^) = (-l)"V27re-x 
A 2n+1 (x) =i(-l) n V2^e^x 2n+1 



4r „2n 

) 

2 



These two expressions are equivalent to (23). 



APPENDIX B: CALCULATION OF THE COEFFICIENTS a AND (3 

We will now calculate coefficients a and j3 of Eq. fl53|). The idea is simply to expand the 
right hand side in the definition of the order parameter, 

/■27T /"OO /"00 

r = / / / cosO -9)p(6,u},n,t)g(n)dndud9. (Bl) 

•/ J —00 J —00 
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as a power series in r. To this end, we fix r and ip in Eq. @, expand its stationary solution 
in powers of r and insert the result in (|Bl"l ). As we do not have a closed formula for the 
stationary solution of (|^), we shall first derive a hierarchy of equations for its coefficients in 
the series in powers of r. Thus we have 

P{9,u,tt;r) = ^ j r, (B2) 



r 



Erf**- (B3) 



n=0 U - 



Then 



rZTV roo roo 

r n = / / cos9G n (9,u,Q)g(fl)dndujde, (B4) 

JO J—oo J—oo 

where we have absorbed ip in the definition of G n . Comparison with eq. (f>3l), establishes 
that ro = T2 = 0, T\ = Ka/2D, and = K 3 j3. From ( |B~2"| ) and (|3|), we obtain 

D d 2 G n 13 <9G n n 



with G-i = 0. Since we are trying to find solutions bifurcating from incoherence, we should 
have G = Po(u,Q,), i.e. the ^-independent incoherent solution @. This directly confirms 
that r = 0. 

The unknowns G n (6, lu, Q) are 27r-periodic in 9, so that we may Fourier expand them as 

oo 

G n (9,cu,n)= Z l n (co,n)e iW , (B6) 

!=— oo 

where 

Z l n (u,Q) = ^ re- iie G n (9,oo,n)d9. (B7) 

Z7T JO 

Here = for I ^ because G is the ^-independent incoherent solution. The condition 
that the probability density function be real yields G n = G n , which in turn implies Z~ l = Z l n . 
Inserting these expressions into ( |B4| ) , we find 

/OO /"OO 
/ Re{Zl)g(Q)dQdu. (B8) 
-OO J — oo 

The unknowns Z l n satisfy the following hierarchy of equations: 

- 2 ^f + ~ [(« - Wl] - * l"Z l n = — — PR ~ Z n-\l (B9) 
m A doo z induo m 1% duo 

The normalization condition for p and the incoherent solution together with ( |B"T| ) imply 

oo 

Z n g(cu)du = 5 n0 . (BIO) 

-oo 
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In order to obtain a, and f3, we should show that r 2 = and calculate t\ and r 3 . We do 
this by means of ( p9|) , with n = 1, 2, 3. 



(a) Case n=l. Eq. (|B9| ) becomes 



D d 2 Z\ Id.. . . - 



m 2 do; 2 m da; 



(Bll) 



where = G is the incoherent solution (|j), and we have used Z 2 = 0. We can solve Eq. 
( |B11| ) by means of the solution of Eq. (|TT|), in which we set A = and replace bi(u, Q) with 
Zl(u,Q), and the right-hand side with i^^dZ^/ 1 duj . Then we obtain 

Z\{^) = -^-e-^--? 

^ !™e^- lV ™°) 2 D p ^dw , s . . 

x V J ~°° 7 («> ), B12 

P =o >»(£+ifi + I>) " 

The previous analogy allows us to calculate ri from ( p8|) and Equations (p0|) , fl22|) and fl23|) 
(with A = 0). The result is 



KD^ mD 



-oo D 2 + ft 2 



mD)P(l + ^ 2 



p=l 



/•oo J 



pi 

1 (B13) 



(b) Case n=2. We now show that r 2 = 0. 

In order to analyze r 2 , it is necessary to consider Z\ y which is the solution to the following 
equation: 

§M + i> - a ^ - ^ - -lllw - ^ < bi4 > 

We shall show that Z\ = Z\ = 0. If this is so, the resulting homogeneous equation can be 
transformed to the parabolic cylinder equation with a quadratic potential having complex 
coefficients. Its solution cannot decay to zero both as ui — > ±oo unless it is identically zero. 
The reason Zf and Z 2 are zero is similar. The equation for Zf is homogeneous (Zq — for 
I ^ 0), and it has a solution C(f2)e~ m<w-f ^ /( 2D ). The normalization condition ( |B10| ) then 
implies C = 0. Z\ again obeys a homogeneous equation which can be transformed into the 
parabolic cylinder equation with a complex potential. The only solution which decays to 
zero as uo — > ±oo is again Z\ = 0. 

(c) Case n=3. 

The equations for Z\ are 
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duj 1 mduj m2i duj 

Dd ^ ■ 1 rf .[(.- fi )z|]-2^zj = -i^ ; , 



2 



rfu; 2 m do; J ~ m 2i du 



D d 2 Z° 2 ,ld ^ 2 K d x . 2K d x 
— 9 - ; — 7T H - = — — Z x - = Im—Z[, 

m z duj z m duj m 1% duo m duj 

D d 2 Z{ Id,. _ . . 1 K d - _ . 

— — f + — (w-fi)Z! \-iuZ] = — Z° (B15) 

m 2 da; 2 mduj LV ; 1J 1 m2icL; V ; 

This system of equations has to be solved numerically in order to obtain the coefficient (3. 
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FIG. 1. Discrete unimoelal frequency distribution: Comparison between analytical (continuous 
line), and numerical (dots) evaluation of the eigenvalue A as a function of m for D = 1 and two 
different values of the coupling strength K. 
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FIG. 2. Discrete unimodal frequency distribution: Eigenvalue A as a function of K for three 
different masses m. Other parameter values and meaning of lines as in Figure 1. Note that curves 
corresponding to different masses all intersect the axis A = at the same point, as expected from 
theory. 
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FIG. 3. Unimodal Lorentzian frequency distribution: Stability diagram of incoherence in the 
parameter space (s,K) for m = 0.2, D = 1. The amplitude of the order parameter as a function 
of time is numerically calculated and displayed at the points marked by: (1) K = 6, e = 1; (2) 
K = 6, e = 1.75; (3) K = 20, e = 4.25; (4) K = 20, e = 4.75. 
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FIG. 5. Discrete bimodal frequency distribution: Stability diagram of incoherence in the pa- 
rameter space (CIq, K) for m = 0.8 and D = 1. The amplitude of the order parameter as a function 
of time is displayed at the points marked by: (1) K = 4.4, f^o = 1-4; (2) K = 5, VLq = 1.4; (3) 
K = 3.6, Q = 0.6; (4) K = 4.4, f\, = 0.6. 
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FIG. 6. Coefficient ci(Oo) corresponding to the approximate mode-coupling stationary solution 
for the discrete bimodal frequency distribution and two different values of m. The solid (m = 0.05) 
and dotted (m = 0.5) lines correspond to the 3-mode approximation, whereas squares and circles 
are obtained from direct numerical simulation. Notice that c\ is even: c\(— $7q) = ci(Oo)- 
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FIG. 7. Coefficient a as a function of the frequency spread e for the stationary solution in 
the case of unimodal Lorentzian frequency distribution. The exact result and that of the 3-mode 
approximation are compared. Other parameter values are D = 1, and (a) m = 0.05; (b) m = 1. 
Note that the approximate result improves as m and e decrease. 
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FIG. 8. Coefficient (3 as a function of the frequency spread e for the stationary solution in 
the case of unimodal Lorentzian frequency distribution. The exact result and that of the 3-mode 
approximation are compared. Other parameter values are D = 1, and (a) m = 0.05; (b) m = 1. 
Note that the approximate result improves as m and e decrease. 
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FIG. 9. Stationary solution, Lorentzian frequency distribution: Relation between frequency 
spread and mass at the critical line f3 = 0. This line separates the regions in parameter space on 
which the synchronization transition from incoherence is either a sub or a supercritical bifurcation. 
The numerically evaluated curve is compared to the approximate result from 3-mode truncation. 
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FIG. 10. Bifurcation diagram for a supercritical bifurcation from incoherence to a synchronized 
stationary solution, corresponding to the unimodal Lorentzian frequency distribution and param- 
eter values of m = 0.05, D = 1, e = 0.5. The amplitude of the order parameter is represented as 
function of K. The continuous line is the analytical approximation and the dots are obtained by 
direct numerical simulation. 
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FIG. 11. Evolution of the amplitude of the order parameter for two different initial conditions, 
unimodal Lorentzian frequency distribution and parameter values m = 0.05, D = 1, e = 5, and 
K = 14.6 (in the region of subcritical bifurcation from incoherence). The simulations illustrate the 
existence of bistability between incoherence and synchronized stationary solution. 
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FIG. 12. Stationary solution, bimodal case: Critical frequency Jig , at which a = 0, as a 
function of mass. Below the critical line, stationary synchronized states bifurcate from incoherence 
at a finite K* . Above the critical line, K* = oo, at which value the branch of synchronized states 
bifurcates subcritically from incoherence. The exact result for the critical line and that obtained 
from 3-mode truncation are also compared. 
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FIG. 13. Stationary solution, bimodal case: Amplitude of the order parameter, r, in the sta- 
tionary synchronized state as a function of the coupling strength, K. We have represented results 
from the 3-mode truncation (lines) and from numerical simulations (dots, squares, circles) for three 
differents values of Oo- Other parameter values are m = 0.1 and D = 1. Notice that K* increases 
with VLq, and it becomes infinity for VLq = 3.5, larger than the critical value S~2g° ps 3.16, as ex- 
pected. This figure illustrates that synchronization issues forth from incoherence as a subcritical 
bifurcation at K* = oo. 



33 




0.15 




(b) 



0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 



FIG. 14. Stationary solution, bimodal case: Coefficient (5 as a function of Qq. We have com- 
pared direct numerical evaluation of (3 (obtained by solving numerically the system of ordinary 
differential equations in Appendix B for the case n=3), and Equation ( |64|) for the discrete bimodal 
frequency distribution. Other parameter values are D = 1 and: (a) m = 0.05, and (b) m = 0.1. 
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FIG. 15. Stationary solution, bimodal case: Critical line at which [3 = in the {Q,Q,m) plane. 
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We have displayed both the direct numerical result and Equation (69). 
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FIG. 16. Stable standing wave solution (SW) for discrete bimodal frequency distribution. We 
have depicted the evolution of the order parameter amplitude r(t) for a large f2o = 15. Other 
parameter values are m = 0.1, D = 1 and K = 5. The inset shows a comparison between the 
numerical solution and the leading-order asymptotic approximation in the high-frequency limit. 
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FIG. 17. Discrete bimodal frequency distribution: Schematic global bifurcation diagram for 
positive values of a and (3 (therefore the bifurcation of the stationary synchronized solution (SS) 
from incoherence is subcritical). Standing and traveling wave solutions are also shown. 
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FIG. 18. Discrete bimodal frequency distribution: Ratio c™ ax /c r Q ax (corresponding to the sta- 
tionary synchronized solution) as a function of n for various mass values and Q = I, D = I. The 
coefficients c n are calculated numerically by integrating the time-dependent equations and waiting 
until a stationary profile is reached. 
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FIG. 19. Unimodal Lorentzian frequency distribution: Comparison between the numerical so- 
lution of a system of N = 20000 Langevin equations, and the numerical method proposed in §V 
containing moments of order 4 or 10. Parameters are m = 0.2, D = 1, e = 1, and K = 8, and we 
have used Q = 15 quadrature nodes. 
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